// Copyright 2011-2020 the Polygon Mesh Processing Library developers.
// Distributed under a MIT-style license, see LICENSE.txt for details.

#include "./SurfaceNormals.h"

namespace zeno {
namespace pmp {

vec3f SurfaceNormals::compute_face_normal(const SurfaceMesh *mesh, int f) {
    int h = mesh->fconn_[f].halfedge_;
    int hend = h;

    auto vpoint = mesh->prim_->attr<vec3f>("pos");

    vec3f p0 = vpoint[mesh->to_vertex(h)];
    h = mesh->next_halfedge(h);
    vec3f p1 = vpoint[mesh->to_vertex(h)];
    h = mesh->next_halfedge(h);
    vec3f p2 = vpoint[mesh->to_vertex(h)];

    if (mesh->next_halfedge(h) == hend) { // face is a triangle
        return normalize(cross(p2 -= p1, p0 -= p1));
    } else { // face is a general polygon
        vec3f n(0, 0, 0);

        // this computes the sum of cross products (area-weighted normals)
        // of triangles generated by inserting the centroid c:
        //   sum_i (p_{i} - c) x (p_{i+1} - c)
        // The point c cancels out, leading to
        //   sum_i (p_{i} x p_{i+1}
        // This vector then has to be normalized.
        for (auto h : mesh->halfedges(f)) {
            n += cross(vpoint[mesh->from_vertex(h)], vpoint[mesh->to_vertex(h)]);
        }

        return normalize(n);
    }
}

vec3f SurfaceNormals::compute_vertex_normal(const SurfaceMesh *mesh, int v) {
    vec3f nn(0.0f);

    if (!mesh->is_isolated(v)) {
        auto &vpoint = mesh->prim_->attr<vec3f>("pos");
        const vec3f p0 = vpoint[v];

        vec3f n;
        vec3f p1, p2;
        float cosine, angle, denom;
        bool is_triangle;

        for (int h : mesh->halfedges(v)) {
            if (mesh->hconn_[h].face_ != PMP_MAX_INDEX) {
                p1 = vpoint[mesh->to_vertex(h)];
                p1 -= p0;
                p2 = vpoint[mesh->from_vertex(mesh->prev_halfedge(h))];
                p2 -= p0;

                // check whether we can robustly compute angle
                denom = sqrt(dot(p1, p1) * dot(p2, p2));
                if (denom > std::numeric_limits<float>::min()) {
                    cosine = dot(p1, p2) / denom;
                    if (cosine < -1.0)
                        cosine = -1.0;
                    else if (cosine > 1.0)
                        cosine = 1.0;
                    angle = acos(cosine);

                    // compute triangle or polygon normal
                    is_triangle = (mesh->next_halfedge(mesh->next_halfedge(mesh->next_halfedge(h))) == h);
                    n = is_triangle ? normalize(cross(p1, p2)) : compute_face_normal(mesh, mesh->hconn_[h].face_);

                    n *= angle;
                    nn += n;
                }
            }
        }

        nn = normalize(nn);
    }

    return nn;
}

void SurfaceNormals::compute_vertex_normals(SurfaceMesh *mesh) {
    auto &vnormal = mesh->prim_->verts.attr<vec3f>("v_normal");
    auto &vdeleted = mesh->prim_->verts.attr<int>("v_deleted");

    for (int v = 0; v < mesh->vertices_size_; ++v) {
        if (mesh->has_garbage_ && vdeleted[v])
            continue;
        vnormal[v] = compute_vertex_normal(mesh, v);
    }
}

} // namespace pmp
} // namespace zeno
